function res = GetWorkingPoint(obj,varargin)
    import MeasurementComponent.Signal.*
    L = log4m.getLogger;
    
    p = inputParser;
    p.addParameter('frequencyZ', 100e3, @(x) x>0);
    p.addParameter('amplitudeZ', 3.0, @(x) x>0 && x<5);
    p.addParameter('frequencyTrans', 20, @(x) x>0);
    p.addParameter('amplitudeTransX', 0.01, @(x) x>0 && x<1);
    p.addParameter('amplitudeTransY', 0.05, @(x) x>0 && x<1);
    p.addParameter('nPoint', 30, @(x) x>=10 && x<=100);
    p.addParameter('nAverage', 10, @(x) x>=1 && x<=100);
    p.addParameter('initialBx0', 0.0, @(x) x>= -3.0 && x <= 3.0);
    p.addParameter('initialBy0', 0.0, @(x) x>= -3.0 && x <= 3.0);
    p.addParameter('quadrant', [0 1], @(x) isnumeric(x) && length(x)==2);
    p.addParameter('transDCscale', 1000, @(x) x>0);
    p.addParameter('transACscale', 100, @(x) x>0);
    p.addParameter('isReport', false, @islogical);
    p.addParameter('isPlot', true, @islogical);
    p.addParameter('description', 'No description.', @ischar);
    p.parse(varargin{:});
    
    %% Prepare environment
    obj.setParametricMagnetometer('frequencyZ', p.Results.frequencyZ, 'amplitudeZ', p.Results.amplitudeZ, 'initialBx0', p.Results.initialBx0, 'initialBy0', p.Results.initialBy0, ...
                                  'bandwidthAC', 5.0*p.Results.frequencyTrans, 'transDCscale', 0, 'transACscale', p.Results.transACscale);
    obj.UnlockTransverseField(); pause(0.5);
    
    %%
    obj.by.setDC(p.Results.initialBy0); 
    obj.bx.reset().setAC('frequency', p.Results.frequencyTrans, 'amplitude', p.Results.amplitudeTransX, 'offset', p.Results.initialBx0);
    bxSig = Signal('bxAmplitude', obj.ziSlave.LockInChannel(4), 'reference', obj.bx);
    obj.bTrans.setupLIAchannel('phase', 0.0+p.Results.quadrant(1)*180); pause(5.0);
    
    bxSig.subscribe_signal(ziDemodTime.R_avg);
    RawPhaseXData = bxSig.getParameterDependence(obj.ziMaster, ziSweepNode.Demod2_PhaseShift, 0.0+p.Results.quadrant(1)*180, 180.0+p.Results.quadrant(1)*180, ...
                                                        'npoint', p.Results.nPoint, 'averaging_sample', p.Results.nAverage);
    RawPhaseXFitting = obj.FitWorkingPoint(RawPhaseXData, 'PhaseRaw'); 
    obj.bTrans.setupLIAchannel('phase', RawPhaseXFitting.phase0-2.0); pause(5.0);

    FinePhaseXData = bxSig.getParameterDependence(obj.ziMaster, ziSweepNode.Demod2_PhaseShift, RawPhaseXFitting.phase0-2.0, RawPhaseXFitting.phase0+2.0, 'npoint', 20, 'averaging_sample', 50);
    FinePhaseXFitting = obj.FitWorkingPoint(FinePhaseXData, 'PhaseFine');
    obj.bx.disable();
    
    %%
    obj.bx.setDC(p.Results.initialBx0); 
    obj.by.reset().setAC('frequency', p.Results.frequencyTrans, 'amplitude', p.Results.amplitudeTransY, 'offset', p.Results.initialBy0,'range',10);
    bySig = Signal('byAmplitude', obj.ziSlave.LockInChannel(4), 'reference', obj.by);
    
    bySig.subscribe_signal(ziDemodTime.R_avg);
    RawPhaseYData = bySig.getParameterDependence(obj.ziMaster, ziSweepNode.Demod2_PhaseShift, 0.0+p.Results.quadrant(2)*180, 180.0+p.Results.quadrant(2)*180, ...
                                                        'npoint', p.Results.nPoint, 'averaging_sample', p.Results.nAverage);
    RawPhaseYFitting = obj.FitWorkingPoint(RawPhaseYData, 'PhaseRaw'); pause(5.0);

    FinePhaseYData = bySig.getParameterDependence(obj.ziMaster, ziSweepNode.Demod2_PhaseShift, RawPhaseYFitting.phase0-2.0, RawPhaseYFitting.phase0+2.0, 'npoint', 20, 'averaging_sample', 50);
    FinePhaseYFitting = obj.FitWorkingPoint(FinePhaseYData, 'PhaseFine');
    obj.by.disable();
    
    %%
    obj.bTrans.setupLIAchannel('phase', FinePhaseYFitting.phase0);
    obj.bTransDCx.setupLIAchannel('phase', FinePhaseYFitting.phase0);
    obj.bTransDCy.setupLIAchannel('phase', FinePhaseXFitting.phase0);
    L.info('NMRGyroscope::GetWorkingPoint', sprintf('PhaseX = %5.3f, PhaseY = %5.3f', FinePhaseXFitting.phase0, FinePhaseYFitting.phase0));
    
    %%
    obj.bx.setDC(p.Results.initialBx0); obj.by.setDC(p.Results.initialBy0);
    transSigGroup = SignalGroup('TransSignals');
    transSigGroup.subscribe_signal(obj.bTransDCx, ziDemodTime.X_avg, [2, 1, 1]);
    transSigGroup.subscribe_signal(obj.bTransDCy, ziDemodTime.X_avg, [2, 1, 2]);
    BiasXData = transSigGroup.getParameterDependence(obj.ziSlave, ziSweepNode.AuxOut1_Offset, p.Results.initialBx0-0.1, p.Results.initialBx0+0.1, 'npoint', 10, 'averaging_sample', 1);
    BiasXFitting = obj.FitWorkingPoint(BiasXData, 'Bias');
    obj.bx.setDC(BiasXFitting.zeropoint);
    L.info('NMRGyroscope::GetWorkingPoint', sprintf('BiasX = %5.3f', BiasXFitting.zeropoint));
    
    %%
    range = [p.Results.initialBy0-0.1, p.Results.initialBy0+0.1];
    obj.by.port.setOutputRange( max(abs(range)) ); pause(0.1);
    transSigGroup.unsubscribe_signal();
    transSigGroup.subscribe_signal(obj.bTransDCy, ziDemodTime.X_avg, [2, 1, 1]);
    transSigGroup.subscribe_signal(obj.bTransDCx, ziDemodTime.X_avg, [2, 1, 2]);
    BiasYData = transSigGroup.getParameterDependence(obj.ziSlave, ziSweepNode.SigOut_Offset, p.Results.initialBy0-0.1, p.Results.initialBy0+0.1, 'npoint', 10, 'averaging_sample', 1);
    BiasYFitting = obj.FitWorkingPoint(BiasYData, 'Bias');
    obj.by.setDC(BiasYFitting.zeropoint);
    L.info('NMRGyroscope::GetWorkingPoint', sprintf('BiasY = %5.3f', BiasYFitting.zeropoint));
    
    %%
    res.timestamp = datetime('now');
    res.signal.rawPhaseX = RawPhaseXData;
    res.signal.rawPhaseY = RawPhaseYData;
    res.signal.finePhaseX = FinePhaseXData;
    res.signal.finePhaseY = FinePhaseYData;
    res.signal.biasX = BiasXData;
    res.signal.biasY = BiasYData;
    res.fitting.rawPhaseX = RawPhaseXFitting;
    res.fitting.rawPhaseY = RawPhaseYFitting;
    res.fitting.finePhaseX = FinePhaseXFitting;
    res.fitting.finePhaseY = FinePhaseYFitting;
    res.fitting.biasX = BiasXFitting;
    res.fitting.biasY = BiasYFitting;
    res.phaseX = FinePhaseXFitting.phase0;
    res.phaseY = FinePhaseYFitting.phase0;
    res.residualX = FinePhaseXFitting.residual;
    res.residualY = FinePhaseYFitting.residual;
    res.biasX = BiasXFitting.zeropoint;
    res.biasY = BiasYFitting.zeropoint;
    res.slopeX = BiasXFitting.slope;
    res.slopeY = BiasYFitting.slope;
    res.exportScaleDC = p.Results.transDCscale;
    res.exportScaleAC = p.Results.transACscale;
    obj.SystemInfo.WorkingPoint = res;
    
    %%
    obj.setParametricMagnetometer('frequencyZ', p.Results.frequencyZ, 'amplitudeZ', p.Results.amplitudeZ, ...
                                  'transDCscale', p.Results.transDCscale, 'transACscale', p.Results.transACscale);
    
    % Plot
    if p.Results.isPlot
        fig = obj.PlotWorkingPoint();
    end
    
    % Record
    if p.Results.isReport
        getReporter;
        rpt.NewRecord(['WorkingPoint_' p.Results.description]);
        if p.Results.isPlot
            rpt.AppendFigure(fig, sprintf('WorkingPoint %s', p.Results.description));
        end
    end
    
    % Store Data
    getSession;
    sess.addData('WorkingPoint', res, p.Results.description);    
end
% 
% function [phase0, fitresult, gof] = swpPhaseFitting(swprResult,varargin)
%     [phaseList, data] = swprResult.getSweeperData();
%     [xData, yData] = prepareCurveData(phaseList, data);
% 
%     amp = max(data); [~, idx] = min(data); phi0 = phaseList(idx);
%     % Set up fittype and options.
%     ft = fittype( 'abs(A*sin(x/180*pi-phi)+B)', 'independent', 'x', 'dependent', 'y' );
%     opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
%     opts.Display = 'Off';
%     opts.Lower =      [amp*0.5 0 0];
%     opts.StartPoint = [amp     0 phi0/180*pi];
%     opts.Upper =      [amp*1.5 0 2*pi];
% 
%     % Fit model to data.
%     [fitresult, gof] = fit( xData, yData, ft, opts );
%     phase0 = fitresult.phi/pi*180;
% end
% 
% function [phase0, residual, fitresult, gof] = swpPhaseFineFitting(swprResult,varargin)
%     [phaseList, data] = swprResult.getSweeperData();
%     [xData, yData] = prepareCurveData(phaseList, data.^2);
%     
%     ft = fittype( 'poly2' );
%     [fitresult, gof] = fit( xData, yData, ft );
%     phase0 = -fitresult.p2/2.0/fitresult.p1;
%     residual = sqrt(fitresult.p3 - fitresult.p2^2/4.0/fitresult.p1);
% end
% 
% function [slope,zeropoint] = swpSlopeFitting(swprResult, varargin)
%     [bTransList, data] = swprResult.getSweeperData();
%     [xData, yData] = prepareCurveData( bTransList, data );
% 
%     ft = fittype( 'poly1' );
%     opts = fitoptions( 'Method', 'LinearLeastSquares' );
% 
%     % Fit model to data.
%     [fitresult, ~] = fit( xData, yData, ft, opts );
%     slope = fitresult.p1;
%     zeropoint = -fitresult.p2./fitresult.p1;
% end
